Regression (Period-based) CLEAN

Intro

Here we provide exploration of LMER models fitted to the 2020 U.S. Congressional election Twitter Dataset (https://doi.org/10.6084/m9.figshare.25523734).

This document was typeset assembled using Quarto (https://quarto.org).

Initialize Workspace

Clear and Init Workspace

Imports

Code
library(ggplot2)
library(GGally)
library(plyr)
library(tidyverse)

Set Directories

Data Load and Inspection

Data Load

Code
ifn0 <- "i0001-politicians-aggregated.rds"
suppressWarnings(rm(list = ls(pattern = "^df")))
df0 <- readr::read_rds(file.path(ifd0, ifn0)) %>%
  dplyr::select(Name, Party, Outcome, Period, Time, Days, Agency) %>%
  dplyr::mutate(
    Period = forcats::fct_recode(
      Period,
      PC = "MI",
      PB = "BE",
      PA = "AR")) %>%
  identity()
df0 %>% dim() %>% cat0()
169997 7 
Code
for (name in names(df0)) {cat0(paste0("    ", name, ", "))}
    Name,  
    Party,  
    Outcome,  
    Period,  
    Time,  
    Days,  
    Agency,  

Skim

Code
set.w(222)
library(skimr)
options(scipen=999)
skimr::skim_format(numeric=list(digits=6))
Warning: Formatting options are now in print() or as function arguments in skim_with().
Code
## digits=4, scientific = FALSE
df0 %>% skimr::skim() %>% dplyr::select(-c(n_missing)) %>% print(numeric=list(digits=2))
── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             169997    
Number of columns          7         
_______________________              
Column type frequency:               
  factor                   4         
  numeric                  3         
________________________             
Group variables            None      

── Variable type: factor ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable complete_rate ordered n_unique top_counts                            
1 Name                      1 FALSE        870 Bar: 361, Ela: 361, Llo: 361, Pra: 361
2 Party                     1 FALSE          2 Dem: 99081, Rep: 70916                
3 Outcome                   1 FALSE          2 win: 115031, los: 54966               
4 Period                    1 FALSE          3 PB: 98660, PA: 47974, PC: 23363       

── Variable type: numeric ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable complete_rate     mean      sd      p0     p25     p50    p75   p100 hist 
1 Time                      1  -0.0680   0.572   -1     -0.55   -0.122  0.428   1    ▇▇▇▆▆
2 Days                      1 -12.2    103.    -180    -99     -22     77     180    ▇▇▇▆▆
3 Agency                    1   0.500    0.271   -1.69   0.343   0.500  0.658   2.37 ▁▁▇▂▁

Check factors order

Code
df2 <- df0
df3 <- df0
summary(df0$Period)
   PC    PB    PA 
23363 98660 47974 
Code
df2$Period <- factor(df2$Period, levels = c("PB", "PC", "PA"))
summary(df2$Period)
   PB    PC    PA 
98660 23363 47974 
Code
df3$Period <- relevel(df0$Period, ref = "PB")
summary(df3$Period)
   PB    PC    PA 
98660 23363 47974 

Prepare for Regression

File names and output directories

## ./data/j1002-model-by-period/i0001_models/fit0x0.extension

Stuff for plotting

Marginalization

Checkups

Save data frame for reference

Code
file <- file.path(dir8, "data-df0.rds")
df0 %>% readr::write_rds(file)
cat0(file)
./data/j1002-model-by-period/i0001_models/data-df0.rds 

Model 01: Null

Fit

Code
model <- "fit01aPd"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (1 | Name) + 1,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

Code
readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)
fit01aPd: [df0] Agency ~ (1 | Name) + 1
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (1 | Name) + 1
   Data: df0
Control: control

REML criterion at convergence: 26631.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-7.9514 -0.5639 -0.0023  0.5683  7.3986 

Random effects:
 Groups   Name        Variance Std.Dev.
 Name     (Intercept) 0.006209 0.0788  
 Residual             0.067519 0.2598  
Number of obs: 169997, groups:  Name, 870

Fixed effects:
              Estimate Std. Error         df t value            Pr(>|t|)    
(Intercept)   0.498618   0.002807 826.426886   177.7 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
cat0(sep2); pp(model); performance::r2(get(model))
--------------------------------------------------------------------- 
fit01aPd: [df0] Agency ~ (1 | Name) + 1
# R2 for Mixed Models

  Conditional R2: 0.084
     Marginal R2: 0.000
Code
cat0(sep2); pp(model); performance::icc(get(model))
--------------------------------------------------------------------- 
fit01aPd: [df0] Agency ~ (1 | Name) + 1
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.084
  Unadjusted ICC: 0.084
Code
cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)
--------------------------------------------------------------------- 
fit01aPd: [df0] Agency ~ (1 | Name) + 1
# ICC by Group

Group |   ICC
-------------
Name  | 0.084

Random Effects

Code
extra <- "9001"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="re")

file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
gg88

Model 02: Time

Fit

Code
model <- "fit02aPd"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

Code
readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)
fit02aPd: [df0] Agency ~ (Time | Name) + Time
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time
   Data: df0
Control: control

REML criterion at convergence: 24941.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.0353 -0.5639 -0.0047  0.5664  7.3362 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Name     (Intercept) 0.006283 0.07927      
          Time        0.004088 0.06394  0.18
 Residual             0.066384 0.25765      
Number of obs: 169997, groups:  Name, 870

Fixed effects:
              Estimate Std. Error         df t value             Pr(>|t|)    
(Intercept)   0.498029   0.002848 816.443563 174.884 < 0.0000000000000002 ***
Time         -0.010230   0.002634 709.116999  -3.884             0.000112 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
     (Intr)
Time 0.187 
Code
cat0(sep2); pp(model); performance::r2(get(model))
--------------------------------------------------------------------- 
fit02aPd: [df0] Agency ~ (Time | Name) + Time
# R2 for Mixed Models

  Conditional R2: 0.102
     Marginal R2: 0.000
Code
cat0(sep2); pp(model); performance::icc(get(model))
--------------------------------------------------------------------- 
fit02aPd: [df0] Agency ~ (Time | Name) + Time
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.102
  Unadjusted ICC: 0.102
Code
cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)
--------------------------------------------------------------------- 
fit02aPd: [df0] Agency ~ (Time | Name) + Time
Model contains random slopes. Cannot compute accurate ICCs by group factors.
# ICC by Group

Group |   ICC
-------------
Name  | 0.085

Model Summary

Code
## pp(model);
tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)
tbl0
Characteristic Beta 95% CI1 p-value
Time -0.01 -0.02, -0.01 <0.001
Name.sd__(Intercept) 0.08

Name.cor__(Intercept).Time 0.18

Name.sd__Time 0.06

Residual.sd__Observation 0.26

1 CI = Confidence Interval

Model Parameters

Code
pp(model);
fit02aPd: [df0] Agency ~ (Time | Name) + Time
Code
parameters::model_parameters(get(model))
# Fixed Effects

Parameter   | Coefficient |       SE |         95% CI | t(169991) |      p
--------------------------------------------------------------------------
(Intercept) |        0.50 | 2.85e-03 | [ 0.49,  0.50] |    174.88 | < .001
Time        |       -0.01 | 2.63e-03 | [-0.02, -0.01] |     -3.88 | < .001

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.08
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |        0.18
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

Plot Model: Estimates

Code
extra <- "0000"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="est") 
gg88 <- gg88 + ggplot2::ylim(-0.1,0.1) + line2 + stuff0
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Code
file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time

Code
terms <- c("Time")
extra <- "0001"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + time0 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time

NULL: Time

Code
extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

pp(model);
fit02aPd: [df0] Agency ~ (Time | Name) + Time
Code
cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.51 | 0.50, 0.51
-0.50 |      0.50 | 0.50, 0.51
 0.00 |      0.50 | 0.49, 0.50
 0.50 |      0.50 | 0.49, 0.50
 1.00 |      0.49 | 0.48, 0.50
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Slope     |       95% CI |     p
--------------------------------
-7.57e-03 | -0.01,  0.00 | 0.004
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.15,)

gg88 <- gg88 + line0 + line1 + time0 + stuff0
## gg88 <- gg88 + coord_cartesian(ylim = c(0.45, 0.55))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

Random Effects

Code
extra <- "9001"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="re")

file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
gg88

Model 03: Time x Period

Fit

Code
model <- "fit03aPd"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time * Period,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

Code
readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time * Period
   Data: df0
Control: control

REML criterion at convergence: 24351.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.1105 -0.5636 -0.0071  0.5640  7.3716 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Name     (Intercept) 0.006226 0.07891      
          Time        0.004116 0.06415  0.18
 Residual             0.066142 0.25718      
Number of obs: 169997, groups:  Name, 870

Fixed effects:
                   Estimate    Std. Error            df t value             Pr(>|t|)    
(Intercept)        0.505972      0.004329   4392.427008 116.881 < 0.0000000000000002 ***
Time              -0.148762      0.016674 151068.334259  -8.922 < 0.0000000000000002 ***
PeriodPB           0.022615      0.003709 169246.914124   6.098        0.00000000108 ***
PeriodPA           0.022052      0.005557 169331.953373   3.968        0.00007248943 ***
Time:PeriodPB      0.193346      0.016754 169082.738995  11.540 < 0.0000000000000002 ***
Time:PeriodPA      0.083065      0.017673 168894.566947   4.700        0.00000260350 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Time   PerdPB PerdPA Tm:PPB
Time        -0.635                            
PeriodPB    -0.699  0.759                     
PeriodPA    -0.461  0.512  0.536              
Time:PerdPB  0.647 -0.976 -0.692 -0.509       
Time:PerdPA  0.616 -0.923 -0.718 -0.759  0.918
Code
cat0(sep2); pp(model); performance::r2(get(model))
--------------------------------------------------------------------- 
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
# R2 for Mixed Models

  Conditional R2: 0.105
     Marginal R2: 0.004
Code
cat0(sep2); pp(model); performance::icc(get(model))
--------------------------------------------------------------------- 
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.101
  Unadjusted ICC: 0.101
Code
cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)
--------------------------------------------------------------------- 
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Model contains random slopes. Cannot compute accurate ICCs by group factors.
# ICC by Group

Group |   ICC
-------------
Name  | 0.085

Model Summary

Code
pp(model);
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Code
tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
Code
tbl0
Characteristic Beta 95% CI1 p-value
Time -0.15 -0.18, -0.12 <0.001
Period


    PB - PC 0.01 0.00, 0.02 0.10
    PA - PC 0.02 0.00, 0.03 0.032
    PA - PB 0.01 0.00, 0.02 0.4
Time * Period


    Time * PB 0.19 0.16, 0.23 <0.001
    Time * PA 0.08 0.05, 0.12 <0.001
Name.sd__(Intercept) 0.08

Name.cor__(Intercept).Time 0.18

Name.sd__Time 0.06

Residual.sd__Observation 0.26

1 CI = Confidence Interval

Model Parameters

Code
pp(model);
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Code
parameters::model_parameters(get(model))
# Fixed Effects

Parameter          | Coefficient |       SE |         95% CI | t(169987) |      p
---------------------------------------------------------------------------------
(Intercept)        |        0.51 | 4.33e-03 | [ 0.50,  0.51] |    116.88 | < .001
Time               |       -0.15 |     0.02 | [-0.18, -0.12] |     -8.92 | < .001
Period [PB]        |        0.02 | 3.71e-03 | [ 0.02,  0.03] |      6.10 | < .001
Period [PA]        |        0.02 | 5.56e-03 | [ 0.01,  0.03] |      3.97 | < .001
Time × Period [PB] |        0.19 |     0.02 | [ 0.16,  0.23] |     11.54 | < .001
Time × Period [PA] |        0.08 |     0.02 | [ 0.05,  0.12] |      4.70 | < .001

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.08
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |        0.18
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

Plot Model: Estimates

Code
pp(model)
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Code
extra <- "0000"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="est")
gg88 <- gg88 + ggplot2::ylim(-0.3,0.3) + line2 + stuff0
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Code
file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time

Code
pp(model);
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Code
extra <- "0001"
terms=c("Time")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + time0 + stuff0 

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Period

Code
pp(model);
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Code
extra <- "0002"
terms=c("Period")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

gg88$data$x <- c(2,1,3)
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Code
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Period

Code
pp(model);
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Code
extra <- "0003"
terms=c("Time", "Period")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time

NULL: Time

Code
pp(model)
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Code
extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.54 | 0.53, 0.55
-0.50 |      0.53 | 0.52, 0.54
 0.00 |      0.53 | 0.52, 0.53
 0.50 |      0.52 | 0.52, 0.53
 1.00 |      0.52 | 0.51, 0.53
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Slope     |       95% CI |     p
--------------------------------
-9.10e-03 | -0.02,  0.00 | 0.025
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05,) 
gg88 <- gg88 + line0 + line1 + time0 + stuff0 
## gg88 <- gg88 + coord_cartesian(ylim = c(0.25, 0.75))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

EF: Period

NULL: Period

Code
pp(model);
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Code
extra <- "1002"
terms <- c("Period")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Period | Predicted |     95% CI
-------------------------------
PC     |      0.52 | 0.51, 0.53
PB     |      0.53 | 0.52, 0.53
PA     |      0.54 | 0.52, 0.55
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
Period | Predicted |     95% CI |      p
----------------------------------------
PB     |      0.53 | 0.52, 0.53 | < .001
PC     |      0.52 | 0.51, 0.53 | < .001
PA     |      0.54 | 0.52, 0.55 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Period

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
## gg88 <- gg88 + line0 + line1 + time0 + stuff0 

gg88 <- gg88 + stuff0 

gg88$data$x <- c(2,1,3)
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Code
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Period

Code
pp(model);
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# Pairwise comparisons

Period |  Contrast |       95% CI |     p
-----------------------------------------
PB-PC  |  9.46e-03 |  0.00,  0.02 | 0.058
PB-PA  | -6.94e-03 | -0.02,  0.00 | 0.173
PC-PA  |     -0.02 | -0.03,  0.00 | 0.036
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Period

NULL: Time x Period

Code
pp(model);
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Code
extra <- "1003"
terms <- c("Time", "Period")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")


cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Period: PC

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.65 | 0.61, 0.69
-0.50 |      0.58 | 0.56, 0.60
 0.00 |      0.51 | 0.50, 0.52
 0.50 |      0.44 | 0.42, 0.45
 1.00 |      0.36 | 0.34, 0.39

Period: PB

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.48 | 0.47, 0.49
-0.50 |      0.51 | 0.50, 0.51
 0.00 |      0.53 | 0.52, 0.54
 0.50 |      0.56 | 0.55, 0.56
 1.00 |      0.58 | 0.57, 0.59

Period: PA

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.59 | 0.57, 0.61
-0.50 |      0.56 | 0.55, 0.58
 0.00 |      0.53 | 0.52, 0.54
 0.50 |      0.50 | 0.49, 0.51
 1.00 |      0.47 | 0.46, 0.48
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Period | Slope |       95% CI |      p
--------------------------------------
PC     | -0.15 | -0.18, -0.12 | < .001
PB     |  0.04 |  0.04,  0.05 | < .001
PA     | -0.07 | -0.08, -0.05 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Period

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + line1 + time0 + stuff0 
gg88 <- gg88 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Time x Period

Code
pp(model);
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# (Average) Linear trend for Time

Period | Contrast |       95% CI |      p
-----------------------------------------
PC-PB  |    -0.19 | -0.23, -0.16 | < .001
PC-PA  |    -0.08 | -0.12, -0.05 | < .001
PB-PA  |     0.11 |  0.10,  0.12 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

Model 04: Time x Period x Outcome

Fit

Code
model <- "fit04aPd"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time * Period * Outcome,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

Code
readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time * Period * Outcome
   Data: df0
Control: control

REML criterion at convergence: 23854.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.1837 -0.5638 -0.0081  0.5631  7.3819 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Name     (Intercept) 0.004807 0.06934       
          Time        0.003341 0.05780  -0.07
 Residual             0.066045 0.25699       
Number of obs: 169997, groups:  Name, 870

Fixed effects:
                                 Estimate    Std. Error            df t value             Pr(>|t|)    
(Intercept)                      0.456679      0.007401  10611.863226  61.707 < 0.0000000000000002 ***
Time                            -0.270178      0.032339 160626.153738  -8.355 < 0.0000000000000002 ***
PeriodPB                         0.058436      0.006870 169298.359477   8.506 < 0.0000000000000002 ***
PeriodPA                        -0.046169      0.011403 168642.624739  -4.049    0.000051518894422 ***
Outcomewinner                    0.077669      0.008954   7522.582122   8.674 < 0.0000000000000002 ***
Time:PeriodPB                    0.327201      0.032499 169002.252327  10.068 < 0.0000000000000002 ***
Time:PeriodPA                    0.250601      0.034945 169220.044400   7.171    0.000000000000746 ***
Time:Outcomewinner               0.157067      0.037712 158776.211635   4.165    0.000031164035167 ***
PeriodPB:Outcomewinner          -0.053725      0.008169 169358.119965  -6.577    0.000000000048304 ***
PeriodPA:Outcomewinner           0.086802      0.013062 169180.936257   6.645    0.000000000030345 ***
Time:PeriodPB:Outcomewinner     -0.177453      0.037926 169245.190328  -4.679    0.000002885388234 ***
Time:PeriodPA:Outcomewinner     -0.213600      0.040515 169254.780591  -5.272    0.000000135027103 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Time   PerdPB PerdPA Otcmwn Tm:PPB Tm:PPA Tm:Otc PrPB:O PrPA:O T:PPB:
Time        -0.715                                                                      
PeriodPB    -0.808  0.769                                                               
PeriodPA    -0.478  0.471  0.514                                                        
Outcomewnnr -0.827  0.591  0.668  0.395                                                 
Time:PerdPB  0.709 -0.985 -0.722 -0.469 -0.586                                          
Time:PerdPA  0.662 -0.912 -0.713 -0.750 -0.547  0.907                                   
Tm:Otcmwnnr  0.613 -0.858 -0.659 -0.404 -0.701  0.845  0.782                            
PrdPB:Otcmw  0.679 -0.646 -0.841 -0.432 -0.770  0.607  0.600  0.765                     
PrdPA:Otcmw  0.417 -0.411 -0.448 -0.873 -0.475  0.410  0.655  0.485  0.519              
Tm:PrdPB:Ot -0.608  0.844  0.619  0.402  0.693 -0.857 -0.778 -0.982 -0.708 -0.482       
Tm:PrdPA:Ot -0.571  0.786  0.615  0.647  0.651 -0.783 -0.863 -0.916 -0.713 -0.753  0.911
Code
cat0(sep2); pp(model); performance::r2(get(model))
--------------------------------------------------------------------- 
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
# R2 for Mixed Models

  Conditional R2: 0.100
     Marginal R2: 0.019
Code
cat0(sep2); pp(model); performance::icc(get(model))
--------------------------------------------------------------------- 
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.083
  Unadjusted ICC: 0.081
Code
cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)
--------------------------------------------------------------------- 
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Model contains random slopes. Cannot compute accurate ICCs by group factors.
# ICC by Group

Group |   ICC
-------------
Name  | 0.067

Model Summary

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
Code
tbl0
Characteristic Beta 95% CI1 p-value
Time -0.27 -0.33, -0.21 <0.001
Period


    PB - PC 0.02 0.00, 0.03 0.007
    PA - PC -0.01 -0.03, 0.01 0.2
    PA - PB -0.03 -0.04, -0.01 <0.001
Outcome


    winner - loser 0.09 0.07, 0.10 <0.001
Time * Period


    Time * PB 0.33 0.26, 0.39 <0.001
    Time * PA 0.25 0.18, 0.32 <0.001
Time * Outcome


    Time * winner 0.16 0.08, 0.23 <0.001
Period * Outcome


    PB * winner -0.05 -0.07, -0.04 <0.001
    PA * winner 0.09 0.06, 0.11 <0.001
Time * Period * Outcome


    Time * PB * winner -0.18 -0.25, -0.10 <0.001
    Time * PA * winner -0.21 -0.29, -0.13 <0.001
Name.sd__(Intercept) 0.07

Name.cor__(Intercept).Time -0.07

Name.sd__Time 0.06

Residual.sd__Observation 0.26

1 CI = Confidence Interval

Model Parameters

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
parameters::model_parameters(get(model))
# Fixed Effects

Parameter                               | Coefficient |       SE |         95% CI | t(169981) |      p
------------------------------------------------------------------------------------------------------
(Intercept)                             |        0.46 | 7.40e-03 | [ 0.44,  0.47] |     61.71 | < .001
Time                                    |       -0.27 |     0.03 | [-0.33, -0.21] |     -8.35 | < .001
Period [PB]                             |        0.06 | 6.87e-03 | [ 0.04,  0.07] |      8.51 | < .001
Period [PA]                             |       -0.05 |     0.01 | [-0.07, -0.02] |     -4.05 | < .001
Outcome [winner]                        |        0.08 | 8.95e-03 | [ 0.06,  0.10] |      8.67 | < .001
Time × Period [PB]                      |        0.33 |     0.03 | [ 0.26,  0.39] |     10.07 | < .001
Time × Period [PA]                      |        0.25 |     0.03 | [ 0.18,  0.32] |      7.17 | < .001
Time × Outcome [winner]                 |        0.16 |     0.04 | [ 0.08,  0.23] |      4.16 | < .001
Period [PB] × Outcome [winner]          |       -0.05 | 8.17e-03 | [-0.07, -0.04] |     -6.58 | < .001
Period [PA] × Outcome [winner]          |        0.09 |     0.01 | [ 0.06,  0.11] |      6.65 | < .001
(Time × Period [PB]) × Outcome [winner] |       -0.18 |     0.04 | [-0.25, -0.10] |     -4.68 | < .001
(Time × Period [PA]) × Outcome [winner] |       -0.21 |     0.04 | [-0.29, -0.13] |     -5.27 | < .001

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.07
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |       -0.07
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

FE: Time

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "0001"
terms=c("Time")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + time0 + stuff0 

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Period

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "0002"
terms=c("Period")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

gg88$data$x <- c(2,1,3)
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Code
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Period

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "0003"
terms=c("Time", "Period")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Outcome

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "0004"
terms=c("Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Outcome

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "0005"
terms=c("Time", "Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + stuff0 + rect3

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Period x Outcome

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "0006"
terms=c("Period", "Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

gg88$data$x <- c(2,1,3) ## CAUTION: THIS WORKS WITH: sjPlot (OR SIMPLER MODEL)
## gg88$data$x <- c(2,2,1,1,3,3) ## CAUTION: THIS WORKS WITH: ggeffects
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Code
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Outcome x Period

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "0007"
terms=c("Time", "Outcome", "Period")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88$data$facet <- factor(
  case_when(
    gg88$data$facet == "PB" ~ "PB",
    gg88$data$facet == "PC" ~ "PC",
    gg88$data$facet == "PA" ~ "PA",
  ),
  levels=c("PB", "PC", "PA"))

gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect0

## gg88 <- gg88 + ggplot2::ylim(0.15,0.85)
## gg88 <- gg88 + ggplot2::coord_cartesian(ylim = c(0.25, 0.75))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time:

NULL: Time

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.54 | 0.53, 0.55
-0.50 |      0.53 | 0.53, 0.54
 0.00 |      0.53 | 0.52, 0.53
 0.50 |      0.52 | 0.52, 0.53
 1.00 |      0.52 | 0.51, 0.53
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Slope |       95% CI |     p
----------------------------
-0.01 | -0.02,  0.00 | 0.005
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05,)

gg88 <- gg88 + line0 + line1 + time0 + stuff0 
## gg88 <- gg88 + coord_cartesian(ylim = c(-0.25, 1.25))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

EF: Period

NULL: Period

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "1002"
terms <- c("Period")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Period | Predicted |     95% CI
-------------------------------
PC     |      0.52 | 0.51, 0.54
PB     |      0.52 | 0.52, 0.53
PA     |      0.52 | 0.51, 0.53
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
Period | Predicted |     95% CI |      p
----------------------------------------
PB     |      0.52 | 0.52, 0.53 | < .001
PC     |      0.52 | 0.51, 0.54 | < .001
PA     |      0.52 | 0.51, 0.53 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Period

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(
  limit_range=FALSE,  
  show_data = FALSE, 
  dot_alpha = 0.05)

gg88 <- gg88 + stuff0
gg88$data$x <- c(2,1,3)
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Code
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Period

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# Pairwise comparisons

Period | Contrast |      95% CI |     p
---------------------------------------
PB-PC  | 3.02e-04 | -0.01, 0.01 | 0.952
PB-PA  | 4.20e-03 | -0.01, 0.01 | 0.874
PC-PA  | 3.90e-03 | -0.01, 0.02 | 0.874
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Period

NULL: Time x Period

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "1003"
terms <- c("Period", "Time")
terms <- c("Time", "Period")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Period: PC

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.67 | 0.63, 0.71
-0.50 |      0.59 | 0.56, 0.61
 0.00 |      0.51 | 0.50, 0.51
 0.50 |      0.43 | 0.41, 0.44
 1.00 |      0.34 | 0.32, 0.37

Period: PB

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.48 | 0.48, 0.49
-0.50 |      0.51 | 0.50, 0.51
 0.00 |      0.53 | 0.52, 0.53
 0.50 |      0.55 | 0.54, 0.56
 1.00 |      0.57 | 0.56, 0.58

Period: PA

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.57 | 0.55, 0.60
-0.50 |      0.55 | 0.53, 0.56
 0.00 |      0.52 | 0.51, 0.53
 0.50 |      0.49 | 0.49, 0.50
 1.00 |      0.46 | 0.46, 0.47
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Period | Slope |       95% CI |      p
--------------------------------------
PC     | -0.27 | -0.33, -0.21 | < .001
PB     |  0.06 |  0.05,  0.07 | < .001
PA     | -0.02 | -0.05,  0.01 | 0.173 
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Period

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
ggsave(file = "file.svg", plot = gg88, width=8, height=6)
gg88

DIFF: Time x Period

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# (Average) Linear trend for Time

Period | Contrast |       95% CI |      p
-----------------------------------------
PC-PB  |    -0.33 | -0.39, -0.26 | < .001
PC-PA  |    -0.25 | -0.32, -0.18 | < .001
PB-PA  |     0.08 |  0.05,  0.11 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Outcome

NULL: Outcome

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "1004"
terms <- c("Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Outcome | Predicted |     95% CI
--------------------------------
loser   |      0.45 | 0.44, 0.46
winner  |      0.52 | 0.51, 0.52
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
Outcome | Predicted |     95% CI |      p
-----------------------------------------
loser   |      0.45 | 0.44, 0.46 | < .001
winner  |      0.52 | 0.51, 0.52 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Outcome

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
ggsave(file = "file.svg", plot = gg88, width=8, height=6)
gg88

DIFF: Outcome

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# Pairwise comparisons

Outcome      | Contrast |       95% CI |      p
-----------------------------------------------
loser-winner |    -0.07 | -0.08, -0.06 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Outcome:

NULL: Time x Outcome

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "1005"
terms <- c("Time", "Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Outcome: loser

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.48 | 0.46, 0.50
-0.50 |      0.48 | 0.47, 0.49
 0.00 |      0.47 | 0.47, 0.48
 0.50 |      0.47 | 0.46, 0.48
 1.00 |      0.47 | 0.45, 0.48

Outcome: winner

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.56 | 0.55, 0.57
-0.50 |      0.55 | 0.54, 0.56
 0.00 |      0.55 | 0.54, 0.55
 0.50 |      0.54 | 0.53, 0.55
 1.00 |      0.53 | 0.52, 0.54
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Outcome | Slope |       95% CI |      p
---------------------------------------
loser   | -0.27 | -0.33, -0.21 | < .001
winner  | -0.11 | -0.15, -0.08 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT Time x Outcome

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect3

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Time x Outcome

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# (Average) Linear trend for Time

Outcome      | Contrast |       95% CI |      p
-----------------------------------------------
loser-winner |    -0.16 | -0.23, -0.08 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Period x Outcome

NULL: Period x Outcome

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "1006"
terms <- c("Period", "Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Outcome: loser

Period | Predicted |     95% CI
-------------------------------
PC     |      0.47 | 0.45, 0.49
PB     |      0.51 | 0.50, 0.52
PA     |      0.41 | 0.39, 0.43

Outcome: winner

Period | Predicted |     95% CI
-------------------------------
PC     |      0.54 | 0.53, 0.55
PB     |      0.53 | 0.53, 0.54
PA     |      0.58 | 0.56, 0.59
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
Period | Outcome | Predicted |     95% CI |      p
--------------------------------------------------
PC     |   loser |      0.47 | 0.45, 0.49 | < .001
PB     |   loser |      0.51 | 0.50, 0.52 | < .001
PA     |   loser |      0.41 | 0.39, 0.43 | < .001
PC     |  winner |      0.54 | 0.53, 0.55 | < .001
PB     |  winner |      0.53 | 0.53, 0.54 | < .001
PA     |  winner |      0.58 | 0.56, 0.59 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Period x Outcome

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)

## gg88$data$x <- c(2,1,3) ## CAUTION: THIS WORKS WITH sjPlot (OR SIMPLER MODEL)
gg88$data$x <- c(2,2,1,1,3,3) ## CAUTION: THIS WORKS WITH: ggeffects
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Code
gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

PLOT: Period x Outcome [low-level]

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "1006"
terms <- c("Period", "Outcome")

type <- "eff"
type <- "pred"

suppressWarnings(rm(list = ls(pattern = "^gg44")))
gg44 <- sjPlot::plot_model(
  get(model), 
  type = type, 
  terms = terms, 
  )

## str(gg44$data)
gg44$data$x <- c(2,2,1,1,3,3)
## gg44 <- gg44 + scale_x_discrete(labels = c("BE", "MI", "ARRR"))
gg44 <- gg44 + scale_x_continuous(
  breaks=c(1,2,3),
  labels=c("PB", "PC", "PA"))
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Code
gg44 <- gg44 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjplot-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg44, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjplot-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg44, width=8, height=6)
gg44

DIFF: Period x Outcome

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# Pairwise comparisons

Period |       Outcome | Contrast |       95% CI |      p
---------------------------------------------------------
PC-PB  |   loser-loser |    -0.04 | -0.05, -0.02 | < .001
PC-PA  |   loser-loser |     0.06 |  0.04,  0.09 | < .001
PC-PC  |  loser-winner |    -0.07 | -0.09, -0.05 | < .001
PC-PB  |  loser-winner |    -0.06 | -0.08, -0.04 | < .001
PC-PA  |  loser-winner |    -0.11 | -0.13, -0.08 | < .001
PB-PA  |   loser-loser |     0.10 |  0.08,  0.12 | < .001
PB-PC  |  loser-winner |    -0.03 | -0.05, -0.02 | < .001
PB-PB  |  loser-winner |    -0.03 | -0.04, -0.01 | < .001
PB-PA  |  loser-winner |    -0.07 | -0.08, -0.05 | < .001
PA-PC  |  loser-winner |    -0.13 | -0.15, -0.11 | < .001
PA-PB  |  loser-winner |    -0.12 | -0.15, -0.10 | < .001
PA-PA  |  loser-winner |    -0.17 | -0.19, -0.14 | < .001
PC-PB  | winner-winner | 5.48e-03 | -0.01,  0.02 | 0.311 
PC-PA  | winner-winner |    -0.04 | -0.05, -0.02 | < .001
PB-PA  | winner-winner |    -0.04 | -0.06, -0.03 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Period x Outcome

NULL: Time x Period x Outcome

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
extra <- "1007"
terms <- c("Time", "Period", "Outcome")
terms <- c("Time", "Outcome", "Period") ## CAUTION
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Outcome: loser
Period: PC

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.72 | 0.65, 0.80
-0.50 |      0.59 | 0.54, 0.63
 0.00 |      0.45 | 0.44, 0.47
 0.50 |      0.32 | 0.30, 0.34
 1.00 |      0.19 | 0.13, 0.24

Outcome: loser
Period: PB

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.45 | 0.44, 0.46
-0.50 |      0.48 | 0.47, 0.49
 0.00 |      0.51 | 0.50, 0.52
 0.50 |      0.54 | 0.53, 0.55
 1.00 |      0.57 | 0.56, 0.59

Outcome: loser
Period: PA

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.43 | 0.38, 0.47
-0.50 |      0.42 | 0.38, 0.45
 0.00 |      0.41 | 0.39, 0.43
 0.50 |      0.40 | 0.39, 0.41
 1.00 |      0.39 | 0.38, 0.41

Outcome: winner
Period: PC

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.64 | 0.60, 0.69
-0.50 |      0.59 | 0.56, 0.61
 0.00 |      0.53 | 0.52, 0.54
 0.50 |      0.48 | 0.46, 0.49
 1.00 |      0.42 | 0.39, 0.45

Outcome: winner
Period: PB

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.50 | 0.49, 0.51
-0.50 |      0.52 | 0.51, 0.52
 0.00 |      0.54 | 0.53, 0.54
 0.50 |      0.56 | 0.55, 0.57
 1.00 |      0.58 | 0.56, 0.59

Outcome: winner
Period: PA

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.65 | 0.62, 0.67
-0.50 |      0.61 | 0.59, 0.63
 0.00 |      0.57 | 0.56, 0.58
 0.50 |      0.54 | 0.53, 0.54
 1.00 |      0.50 | 0.49, 0.51
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Outcome | Period | Slope |       95% CI |      p
------------------------------------------------
loser   |     PC | -0.27 | -0.33, -0.21 | < .001
loser   |     PB |  0.06 |  0.05,  0.07 | < .001
loser   |     PA | -0.02 | -0.05,  0.01 | 0.173 
winner  |     PC | -0.11 | -0.15, -0.08 | < .001
winner  |     PB |  0.04 |  0.03,  0.05 | < .001
winner  |     PA | -0.08 | -0.09, -0.06 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Period x Outcome

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)

gg88$data$facet <- factor(
  case_when(
    gg88$data$facet == "PB" ~ "PB",
    gg88$data$facet == "PC" ~ "PC",
    gg88$data$facet == "PA" ~ "PA",
  ),
  levels=c("PB", "PC", "PA"))

gg88 <- gg88 + line0 + line1 + time0 + stuff0 

gg88 <- gg88 + rect0

## gg88 <- gg88 + ggplot2::ylim(0.15,0.85)

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=12, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=12, height=6)
gg88

PLOT: Time x Period x Outcome (low-level)

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
          get(model),
          type="emm",
          terms=terms) 

gg88$data$facet <- factor(
  case_when(
    gg88$data$facet == "PB" ~ "PB",
    gg88$data$facet == "PC" ~ "PC",
    gg88$data$facet == "PA" ~ "PA",
  ),
  levels=c("PB", "PC", "PA"))

gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect0
## gg88 <- gg88 + ggplot2::ylim(0.15,0.85)

file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjPlot-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=12, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjPlot-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=12, height=6)
gg88

DIFF: Time x Period x Outcome

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# (Average) Linear trend for Time

Outcome       | Period | Contrast |       95% CI |      p
---------------------------------------------------------
loser-loser   |  PC-PB |    -0.33 | -0.39, -0.26 | < .001
loser-loser   |  PC-PA |    -0.25 | -0.32, -0.18 | < .001
loser-winner  |  PC-PC |    -0.16 | -0.23, -0.08 | < .001
loser-winner  |  PC-PB |    -0.31 | -0.37, -0.24 | < .001
loser-winner  |  PC-PA |    -0.19 | -0.26, -0.13 | < .001
loser-loser   |  PB-PA |     0.08 |  0.05,  0.11 | < .001
loser-winner  |  PB-PC |     0.17 |  0.13,  0.21 | < .001
loser-winner  |  PB-PB |     0.02 |  0.01,  0.03 | 0.005 
loser-winner  |  PB-PA |     0.13 |  0.11,  0.15 | < .001
loser-winner  |  PA-PC |     0.09 |  0.05,  0.14 | < .001
loser-winner  |  PA-PB |    -0.06 | -0.09, -0.03 | < .001
loser-winner  |  PA-PA |     0.06 |  0.02,  0.09 | < .001
winner-winner |  PC-PB |    -0.15 | -0.19, -0.11 | < .001
winner-winner |  PC-PA |    -0.04 | -0.08,  0.00 | 0.071 
winner-winner |  PB-PA |     0.11 |  0.10,  0.13 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

CHECK: Model Assumptions

Code
pp(model);
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome
Code
performance::check_model(get(model))

Code
ggsave(
  file = "performance-test.png", 
  ## file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  ## plot = gg88,
  ## plot = last_plot(),
  width=8,
  height=12)
Warning: Failed to fit group -1.
Failed to fit group -1.
Caused by error in `predLoess()`:
! workspace required (43350042532) is too large probably because of setting 'se = TRUE'.
Code
ggsave(
  file = "performance-test.png", 
  ## file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  ## plot = gg88,
  ## plot = last_plot(),
  width=8,
  height=12)
Warning: Failed to fit group -1.
Failed to fit group -1.
Caused by error in `predLoess()`:
! workspace required (43350042532) is too large probably because of setting 'se = TRUE'.
Code
ls()
 [1] "cat0"                     "check_pkgs"               "control"                  "create_factor_mappings"   "df0"                      "df2"                      "df3"                      "dir8"                    
 [9] "extra"                    "fbase"                    "file"                     "fit01aPd"                 "fit02aPd"                 "fit03aPd"                 "fit04aPd"                 "fpath"                   
[17] "generate_permutations"    "get_stars"                "gg44"                     "gg88"                     "ggpairs_lower_fun"        "hasConverged"             "ifd0"                     "ifn0"                    
[25] "jiko_count_na"            "line0"                    "line1"                    "line2"                    "model"                    "name"                     "obj_has_name"             "ofd0"                    
[33] "ofd2"                     "pp"                       "pred0"                    "pred0tab"                 "rect0"                    "rect2"                    "rect3"                    "REML"                    
[41] "replace_factor_levels"    "replace_outliers_with_NA" "rescale01"                "rxy0"                     "rxy2"                     "rxy3"                     "sep0"                     "sep1"                    
[49] "sep2"                     "set.w"                    "stuff0"                   "tbl0"                     "terms"                    "test0"                    "test0tab"                 "test2"                   
[57] "test2tab"                 "time0"                    "type"                     "var0"                     "write_model_info"        

Auxiliary Checkups

Variables influence / importance in fit04aPd

Code
## TODO CHANGE TITLE
formula(fit04aPd)
Agency ~ (Time | Name) + Time * Period * Outcome
Code
eff_size0 <- effectsize::eta_squared(fit04aPd)
eff_size2 <- eff_size0 %>%
  dplyr::mutate(Interpret = effectsize::interpret_eta_squared(Eta2_partial))

performance::print_html(eff_size2)
Effect Size for ANOVA
Parameter Eta2 (partial) 95% CI Interpret
Time 2.12e-03 [0.00, 1.00] very small
Period 5.13e-04 [0.00, 1.00] very small
Outcome 0.07 [0.05, 1.00] medium
Time:Period 1.65e-03 [0.00, 1.00] very small
Time:Outcome 9.19e-05 [0.00, 1.00] very small
Period:Outcome 1.08e-03 [0.00, 1.00] very small
Time:Period:Outcome 1.71e-04 [0.00, 1.00] very small

List models (for Model Hierarchy / Hierarchies)

Code
temp <- lapply(ls(pattern="^fit\\d+"), pp)
fit01aPd: [df0] Agency ~ (1 | Name) + 1
fit02aPd: [df0] Agency ~ (Time | Name) + Time
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

Compare models using performance package with score

Code
perf0 <- performance::compare_performance(
  fit01aPd, # [df0] Agency ~ (1 | Name) + 1
  fit02aPd, # [df0] Agency ~ (Time | Name) + Time
  fit03aPd, # [df0] Agency ~ (Time | Name) + Time * Period
  fit04aPd, # [df0] Agency ~ (Time | Name) + Time * Period * Outcome
  ## CAUTION: COMMA
  rank = TRUE, verbose = FALSE)
perf0 %>% performance::print_html()
Comparison of Model Performance Indices
Name Model R2 (cond.) R2 (marg.) ICC RMSE Sigma AIC weights AICc weights BIC weights Performance-Score
fit04aPd lmerModLmerTest 0.10 0.02 0.08 0.26 0.26 1.00 1.00 1.00 84.36%
fit03aPd lmerModLmerTest 0.11 4.32e-03 0.10 0.26 0.26 4.74e-116 4.74e-116 5.77e-103 51.23%
fit02aPd lmerModLmerTest 0.10 4.62e-04 0.10 0.26 0.26 5.28e-250 5.29e-250 3.40e-228 43.11%
fit01aPd lmerModLmerTest 0.08 0.00 0.08 0.26 0.26 0.00e+00 0.00e+00 0.00e+00 1.03%
NA
Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- perf0 %>% plot()
ggsave(
  file = file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  plot = gg88,
  width=8,
  height=6)

gg88

Compare less models

Code
perf2 <- performance::compare_performance(
  ## fit01aPd, # [df0] Agency ~ (1 | Name) + 1
  ## fit02aPd, # [df0] Agency ~ (Time | Name) + Time
  fit03aPd, # [df0] Agency ~ (Time | Name) + Time * Period
  fit04aPd, # [df0] Agency ~ (Time | Name) + Time * Period * Outcome
  ## CAUTION: COMMA
  rank = TRUE, verbose = FALSE)
perf2 %>% performance::print_html()
Comparison of Model Performance Indices
Name Model R2 (cond.) R2 (marg.) ICC RMSE Sigma AIC weights AICc weights BIC weights Performance-Score
fit04aPd lmerModLmerTest 0.10 0.02 0.08 0.26 0.26 1.00 1.00 1.00 75.00%
fit03aPd lmerModLmerTest 0.11 4.32e-03 0.10 0.26 0.26 4.74e-116 4.74e-116 5.77e-103 25.00%
NA
Code
suppressWarnings(rm(list = ls(pattern = "^gg44")))
gg44 <- perf2 %>% plot()
ggsave(
  file = file.path(dir8, "summary-combined-models.perf2.compare-03-score.png"),
  plot = gg44,
  width=8,
  height=6)

gg44  

Comment

Add interpretation of the difference between best fit to data and best predictive power and overfitting (aspecially in the context of the last model).

Code
gg44

Performance Table Sorted by R2_conditional

Code
perf0 %>% dplyr::arrange(desc(R2_conditional)) %>% performance::print_html()
Comparison of Model Performance Indices
Name Model R2 (cond.) R2 (marg.) ICC RMSE Sigma AIC weights AICc weights BIC weights Performance-Score
fit03aPd lmerModLmerTest 0.11 4.32e-03 0.10 0.26 0.26 4.74e-116 4.74e-116 5.77e-103 51.23%
fit02aPd lmerModLmerTest 0.10 4.62e-04 0.10 0.26 0.26 5.28e-250 5.29e-250 3.40e-228 43.11%
fit04aPd lmerModLmerTest 0.10 0.02 0.08 0.26 0.26 1.00 1.00 1.00 84.36%
fit01aPd lmerModLmerTest 0.08 0.00 0.08 0.26 0.26 0.00e+00 0.00e+00 0.00e+00 1.03%
NA

Performance Table Sorted by R2_marginal

Code
perf0 %>% dplyr::arrange(desc(R2_marginal)) %>% performance::print_html()
Comparison of Model Performance Indices
Name Model R2 (cond.) R2 (marg.) ICC RMSE Sigma AIC weights AICc weights BIC weights Performance-Score
fit04aPd lmerModLmerTest 0.10 0.02 0.08 0.26 0.26 1.00 1.00 1.00 84.36%
fit03aPd lmerModLmerTest 0.11 4.32e-03 0.10 0.26 0.26 4.74e-116 4.74e-116 5.77e-103 51.23%
fit02aPd lmerModLmerTest 0.10 4.62e-04 0.10 0.26 0.26 5.28e-250 5.29e-250 3.40e-228 43.11%
fit01aPd lmerModLmerTest 0.08 0.00 0.08 0.26 0.26 0.00e+00 0.00e+00 0.00e+00 1.03%
NA

Interpret R2

Code
model <- "fit01aPd" # [df0] Agency ~ (1 | Name) + 1
model <- "fit02aPd" # [df0] Agency ~ (Time | Name) + Time
model <- "fit03aPd" # [df0] Agency ~ (Time | Name) + Time * Period
model <- "fit04aPd" # [df0] Agency ~ (Time | Name) + Time * Period * Outcome

cat0(effectsize::interpret_r2(performance::r2(get(model))$R2_conditional, rules="cohen1988"))
weak 
Code
cat0(effectsize::interpret_r2(performance::r2(get(model))$R2_marginal, rules="cohen1988"))
very weak 

Plot models: Part 1: Basic Models

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_models(
  ## CAUTION the null model can not be used here 
  ## Thus to keep the numbers consistent I have 
  ## used model 02 as an input twice
  fit02aPd, # [df0] Agency ~ (Time | Name) + Time 
  fit02aPd, # [df0] Agency ~ (Time | Name) + Time
  fit03aPd, # [df0] Agency ~ (Time | Name) + Time * Period
  fit04aPd, # [df0] Agency ~ (Time | Name) + Time * Period * Outcome
  spacing=1, dot.size=1
)

ggsave(
  file = file.path(dir8, "summary-all-models-part-01.png"),
  plot = gg88,
  width=5,
  height=4)

gg88

Tab Models

The above `plot_models` seems to be a lot more illustrative

Code
# , eval = TRUE, results='hide'
library(sjPlot)
library(sjmisc)

Attaching package: 'sjmisc'
The following object is masked from 'package:skimr':

    to_long
The following object is masked from 'package:purrr':

    is_empty
The following object is masked from 'package:tidyr':

    replace_na
The following object is masked from 'package:tibble':

    add_case
Code
library(sjlabelled)

Attaching package: 'sjlabelled'
The following object is masked from 'package:forcats':

    as_factor
The following object is masked from 'package:dplyr':

    as_label
The following object is masked from 'package:ggplot2':

    as_label
Code
tab0 <- sjPlot::tab_model(
  fit01aPd, # [df0] Agency ~ (1 | Name) + 1
  fit02aPd, # [df0] Agency ~ (Time | Name) + Time
  fit03aPd, # [df0] Agency ~ (Time | Name) + Time * Period
  fit04aPd, # [df0] Agency ~ (Time | Name) + Time * Period * Outcome
    show.reflvl = TRUE,
    show.intercept = TRUE,
    p.style = "numeric_stars",
    file = file.path(dir8, "summary-all-models-part-01-SEE-PNG-images-for-a-better-overview.html"))
tab0
  Agency Agency Agency Agency
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.50 *** 0.49 – 0.50 <0.001 0.50 *** 0.49 – 0.50 <0.001 0.51 *** 0.50 – 0.51 <0.001 0.46 *** 0.44 – 0.47 <0.001
PC Reference Reference Reference Reference
PB 0.02 *** 0.02 – 0.03 <0.001 0.06 *** 0.04 – 0.07 <0.001
PA 0.02 *** 0.01 – 0.03 <0.001 -0.05 *** -0.07 – -0.02 <0.001
PeriodPA:Outcomewinner 0.09 *** 0.06 – 0.11 <0.001
loser Reference Reference Reference Reference
PeriodPB:Outcomewinner -0.05 *** -0.07 – -0.04 <0.001
winner 0.08 *** 0.06 – 0.10 <0.001
Time -0.01 *** -0.02 – -0.01 <0.001 -0.15 *** -0.18 – -0.12 <0.001 -0.27 *** -0.33 – -0.21 <0.001
Time:Outcomewinner 0.16 *** 0.08 – 0.23 <0.001
Time:PeriodPA 0.08 *** 0.05 – 0.12 <0.001 0.25 *** 0.18 – 0.32 <0.001
Time:PeriodPA:Outcomewinner -0.21 *** -0.29 – -0.13 <0.001
Time:PeriodPB 0.19 *** 0.16 – 0.23 <0.001 0.33 *** 0.26 – 0.39 <0.001
Time:PeriodPB:Outcomewinner -0.18 *** -0.25 – -0.10 <0.001
Random Effects
σ2 0.07 0.07 0.07 0.07
τ00 0.01 Name 0.01 Name 0.01 Name 0.00 Name
τ11   0.00 Name.Time 0.00 Name.Time 0.00 Name.Time
ρ01   0.18 Name 0.18 Name -0.07 Name
ICC 0.08 0.10 0.10 0.08
N 870 Name 870 Name 870 Name 870 Name
Observations 169997 169997 169997 169997
Marginal R2 / Conditional R2 0.000 / 0.084 0.000 / 0.102 0.004 / 0.105 0.019 / 0.100
* p<0.05   ** p<0.01   *** p<0.001
Code
## knitr::html
## papaja::
## rmarkdown::render
## rmarkdown::render(tab0$knitr)
## tab0$knitr

Tabulate models (publication ready)

Code
tab0 <- sjPlot::tab_model(
  fit01aPd, # [df0] Agency ~ (1 | Name) + 1
  fit02aPd, # [df0] Agency ~ (Time | Name) + Time
  fit03aPd, # [df0] Agency ~ (Time | Name) + Time * Period
  fit04aPd, # [df0] Agency ~ (Time | Name) + Time * Period * Outcome
    show.reflvl = FALSE,
    show.intercept = TRUE,
    p.style = "numeric_stars",
    file = file.path(dir8, "summary-all-models-part-02-publication.html"))
tab0
  Agency Agency Agency Agency
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.50 *** 0.49 – 0.50 <0.001 0.50 *** 0.49 – 0.50 <0.001 0.51 *** 0.50 – 0.51 <0.001 0.46 *** 0.44 – 0.47 <0.001
Time -0.01 *** -0.02 – -0.01 <0.001 -0.15 *** -0.18 – -0.12 <0.001 -0.27 *** -0.33 – -0.21 <0.001
Period [PB] 0.02 *** 0.02 – 0.03 <0.001 0.06 *** 0.04 – 0.07 <0.001
Period [PA] 0.02 *** 0.01 – 0.03 <0.001 -0.05 *** -0.07 – -0.02 <0.001
Time × Period [PB] 0.19 *** 0.16 – 0.23 <0.001 0.33 *** 0.26 – 0.39 <0.001
Time × Period [PA] 0.08 *** 0.05 – 0.12 <0.001 0.25 *** 0.18 – 0.32 <0.001
Outcome [winner] 0.08 *** 0.06 – 0.10 <0.001
Time × Outcome [winner] 0.16 *** 0.08 – 0.23 <0.001
Period [PB] × Outcome
[winner]
-0.05 *** -0.07 – -0.04 <0.001
Period [PA] × Outcome
[winner]
0.09 *** 0.06 – 0.11 <0.001
(Time × Period [PB]) ×
Outcome [winner]
-0.18 *** -0.25 – -0.10 <0.001
(Time × Period [PA]) ×
Outcome [winner]
-0.21 *** -0.29 – -0.13 <0.001
Random Effects
σ2 0.07 0.07 0.07 0.07
τ00 0.01 Name 0.01 Name 0.01 Name 0.00 Name
τ11   0.00 Name.Time 0.00 Name.Time 0.00 Name.Time
ρ01   0.18 Name 0.18 Name -0.07 Name
ICC 0.08 0.10 0.10 0.08
N 870 Name 870 Name 870 Name 870 Name
Observations 169997 169997 169997 169997
Marginal R2 / Conditional R2 0.000 / 0.084 0.000 / 0.102 0.004 / 0.105 0.019 / 0.100
* p<0.05   ** p<0.01   *** p<0.001
Code
# save.image(file=file.path(dir8, "session.RData"))
# load(file.path(dir8, "session.RData")

Report

Code
result <- report::report(fit04aPd)
Code
print(result)
We fitted a linear mixed model (estimated using REML and Nelder-Mead optimizer) to predict Agency with Time, Period and Outcome (formula: Agency ~ Time * Period * Outcome). The model included Time as random effects
(formula: ~Time | Name). The model's total explanatory power is weak (conditional R2 = 0.10) and the part related to the fixed effects alone (marginal R2) is of 0.02. The model's intercept, corresponding to Time = 0,
Period = PC and Outcome = loser, is at 0.46 (95% CI [0.44, 0.47], t(169981) = 61.71, p < .001). Within this model:

  - The effect of Time is statistically significant and negative (beta = -0.27, 95% CI [-0.33, -0.21], t(169981) = -8.35, p < .001; Std. beta = -0.57, 95% CI [-0.70, -0.44])
  - The effect of Period [PB] is statistically significant and positive (beta = 0.06, 95% CI [0.04, 0.07], t(169981) = 8.51, p < .001; Std. beta = 0.13, 95% CI [0.07, 0.20])
  - The effect of Period [PA] is statistically significant and negative (beta = -0.05, 95% CI [-0.07, -0.02], t(169981) = -4.05, p < .001; Std. beta = -0.23, 95% CI [-0.33, -0.14])
  - The effect of Outcome [winner] is statistically significant and positive (beta = 0.08, 95% CI [0.06, 0.10], t(169981) = 8.67, p < .001; Std. beta = 0.25, 95% CI [0.17, 0.33])
  - The effect of Time × Period [PB] is statistically significant and positive (beta = 0.33, 95% CI [0.26, 0.39], t(169981) = 10.07, p < .001; Std. beta = 0.69, 95% CI [0.56, 0.83])
  - The effect of Time × Period [PA] is statistically significant and positive (beta = 0.25, 95% CI [0.18, 0.32], t(169981) = 7.17, p < .001; Std. beta = 0.53, 95% CI [0.38, 0.67])
  - The effect of Time × Outcome [winner] is statistically significant and positive (beta = 0.16, 95% CI [0.08, 0.23], t(169981) = 4.16, p < .001; Std. beta = 0.33, 95% CI [0.18, 0.49])
  - The effect of Period [PB] × Outcome [winner] is statistically significant and negative (beta = -0.05, 95% CI [-0.07, -0.04], t(169981) = -6.58, p < .001; Std. beta = -0.15, 95% CI [-0.23, -0.08])
  - The effect of Period [PA] × Outcome [winner] is statistically significant and positive (beta = 0.09, 95% CI [0.06, 0.11], t(169981) = 6.65, p < .001; Std. beta = 0.37, 95% CI [0.26, 0.48])
  - The effect of (Time × Period [PB]) × Outcome [winner] is statistically significant and negative (beta = -0.18, 95% CI [-0.25, -0.10], t(169981) = -4.68, p < .001; Std. beta = -0.37, 95% CI [-0.53, -0.22])
  - The effect of (Time × Period [PA]) × Outcome [winner] is statistically significant and negative (beta = -0.21, 95% CI [-0.29, -0.13], t(169981) = -5.27, p < .001; Std. beta = -0.45, 95% CI [-0.62, -0.28])

Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using a Wald t-distribution approximation.